﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using GraphEntities;
using NumericalData.Interpolation.Data;

namespace NumericalData.Interpolation.Data.Methods
{
    abstract public class InterpolatedFunction
    {
        protected ExpressionContainer _container;
        public ExpressionContainer ExpressionContainer { get; set; }

        public InterpolatedFunction()
        {
        }

        public InterpolatedFunction(ExpressionContainer exp)
        {
            _container = exp.IsValid ? exp : null;
        }

        public Points CalculateValues(double[] xValues)
        {
            var points = new Points();
            foreach (double x in xValues)
            {
                points.Add(new Point(x,Invoke(x)));
            }
            return points;
        }


        public abstract double Invoke(double x);

    }

    public class StirlingFunction : InterpolatedFunction
    {
        private double _h;


        public StirlingFunction(ExpressionContainer exp)
            : base(exp)
        {
            _h = exp.Grid.Step.Value;
        }

        public override double Invoke(double x)
        {
            var order = _container.MethodOrder;
            var points = _container.Grid.Values;
            var x0Index = CalculateClosestPoint(x);
            var x0 = points[x0Index].X;
            var y0 = points[x0Index].Y;
            if (x0Index - order < 0 || x0Index + order > points.Count - 1)
                throw new Exception("Can't calculate interpolation function with this order in this point");
            var q = (x - x0) / _h;
            double acc = y0;
            var deltas = new DeltasContainer(_container.Grid, x0Index);
            for (var i = 1; i <= order; ++i)
            {
                Console.WriteLine(String.Format("iteration {0}", i));
                var ei = q * Ei(q, i);
                var ai = 2 * i * ei * (deltas.GetDelta(-i, 2 * i - 1) + deltas.GetDelta(-i + 1, 2 * i - 1)) / 2;
                var bi = ei * deltas.GetDelta(-i, 2 * i);
                acc += ai + bi;
            }
            return acc;
        }

        private int CalculateClosestPoint(double x)
        {
            if (x > _container.IntervalStop + Double.Epsilon ||
                x < _container.IntervalStart - Double.Epsilon)
                throw new Exception("CalculateClosestPoint: x is out of range");
            var index = 0;
            var measure = Math.Abs(x - _container.Grid.Values[index].X);

            while (measure > Math.Abs(x - _container.Grid.Values[++index].X))
                measure = Math.Abs(x - _container.Grid.Values[index].X);

            return index - 1;
        }

        private double Ei(double q, int i)
        {
            var acc = q / 2;
            for (var k = 1; k <= i - 1; k++)
            {
                acc *= (q * q - k * k) / ((1 + 2 * k) * (2 + 2 * k));
            }
            return acc;
        }
    }
}